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Abstract. Criterion for contacting is critically important for the Generalized 
Interpolation Material Point (GIMP) method. We present an improved criterion 
by adding a switching function. With the method dynamical response of high 
melting explosive (HMX) with cavities under shock is investigated. The physical 
model used in the present work is an clastic-to-plastic and thermal-dynamical model 
with Mie-Griineissen equation of state. We mainly concern the influence of various 
parameters, including the impacting velocity v, cavity size R, etc, to the dynamical 
and thermodynamical behaviors of the material. For the colliding of two bodies with a 
cavity in each, a secondary impacting is observed. Correspondingly, the separation 
distance D of the two bodies has a maximum value D max in between the initial 
and second impacts. When the initial impacting velocity v is not large enough, the 
cavity collapses in a nearly symmetric fashion, the maximum separation distance D max 
increases with v. When the initial shock wave is strong enough to collapse the cavity 
asymmetrically along the shock direction, the variation of D max with v does not show 
monotonic behavior. Our numerical results show clear indication that the existence of 
cavities in explosive helps the creation of "hot spots" . 
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1. Introduction 

Cavitation phenomena are ubiquitous in nature, ranging from solid, liquid to plasma. 
Cavity creation and collapse play a very important role in a large amount of industrial 
processes, like erosion of materials, ignition of explosive, comminution of kidney stones, 
etc [1-3]. The cavitation phenomena may occur in a mesoscopic scale, but generally 
have a profound influence on the response in the macroscopic scale. As for explosives, 
the cavity collapse may result in jetting phenomena and consequently lead to a much 
higher temperature increase of the "hot spots". Therefore, it has a potential to start 
local reaction leading to partial decomposition or run to full detonation [2]. On the one 
hand, cavitation acts as a means of increasing the sensitivity of explosive for ignition, but 
on the other hand, it represents a potential safety problem for handling such materials. 
The interest in both considerations inspired an amount of theoretical, experimental and 
numerical work aiming to study the influence of cavities to the explosives [2,4-14]. 

During the collapse procedure explosive material under consideration is highly 
distorted and the boundaries should be tracked, which proposes a very high requirement 
to the numerical methods. In the early studies, conventional Eulerian and/or Lagrangian 
schemes were used [5-9]. When conventional Lagrangian methods are used, the mesh 
distortion associated with large deformations reduces the accuracy and will probably 
terminate the computations; ultimately, a re-meshing treatment is required, which is 
not a straightforward task and generally inefficient. It is known that the conventional 
Eulerian code is not convenient to track the boundaries. To continue, several mixed 
methods have been proposed to combine the merits of the two methods and overcome 
their drawbacks. Among them, the Arbitrary Lagrange- Euler (ALE) method [15] is 
a typical one. Several researchers have used the ALE method to study the collapse 
of cavities in explosive materials [16]. The particle-in-cell (PIC) is a second typical 
mixed method [17]. It also contributed to the research work on inhomogeneous plastic- 
bonded-explosives [18]. Additionally, a number of "meshless methods" have also been 
developed. To overcome the difficulty due to the distortion of mesh, the meshless method 
uses a series of discretized lagrangian points to construct the shape function. Some 
meshless methods, such as the Smooth Particle Hydrodynamics(SPH) [19-21], Dual 
Particle Dynamics(DPD) [22], etc, have been applied in the research of cavity collapse 
in explosive materials [22-24]. It should be mentioned that the microscopic Molecular 
Dynamics(MD) has also brought some new understanding of the physics and chemistry 
involved [12-14]. Yet, even with the most powerful computers in nowadays, the MD 
simulation is still far from reaching the practical spatial and temporal scales of real 
experiments [25]. 

The physical models used in previous numerical studies are mainly fluidic ones. 
Such models ignore many important characteristics of solid, like plastic strain, 
hardening, and effects relevant to the deformation history. In this study we revisit the 
Material Point Method(MPM) [26], recently developed in the field of computational solid 
physics, by presenting a new contact criterion, and then use it to study the dynamical 
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response of high melting explosive (HMX) with cavities under shock. 

The MPM is a descendant of the PIC extended to solid mechanics [27-31]. 
Compared with the above-mentioned methods, the MPM provides a robust and efficient 
treatment of large deformation issues and a convenient framework for modeling contact 
between large numbers of contacting bodies. In MPM, each body is discretized by a 
collection of Lagrangian material points carrying all information required to advance the 
solution; the Eulerian background computational mesh is used to solve the governing 
equations and the mesh solution is then used to update the information on material 
points. Specifically, at each time step, the MPM calculations consist of two parts: a 
Lagrangian part and a convective one. Firstly, the computational mesh deforms with 
the body, and is used to determine the strain increment and the stresses carried by the 
particles. Then, the new position of the computational mesh is chosen (particularly, it 
may be the previous one); and the velocity field is mapped from the material points to 
the mesh nodes. This method has been applied to handle engineering problems with 
large strain [32,33] and/or dynamical energy release rate [34], fracture in heterogeneous 
material [35], dynamics failure [36,37], hypervelocity impact [38], thin membranes [39], 
granular materials [18,40-42], etc. The Generalized Interpolation Material Point(GIMP) 
method uses a variational form and a Petrov-Galerkin discretization scheme to overcome 
the numerical noise of previous MPM [26]. 

We study not only the collapse of cavities but also the deformation and 
dynamical response of cavities before collapsing under shock. These are important 
for understanding better the ignition of explosive. The elastic-to-plastic and thermal- 
dynamics model with the Mie-Griineissen equation of state consulting the Rankine- 
Hugoniot curve is used to simulate the mechanical and thermal behaviors of the material. 
The model used here is more suitable for simulating the dynamical response of the solid 
material under shock than previous ones. 

This paper is organized as follows. In Sect. 2 we briefly describe the generalized 
interpolation material point method where the contact algorithm is improved. In Sect. 
3 the new scheme is validated by several benchmark tests and then is used to study 
cavitated high melting explosive. Sect. 4 concludes the present paper. 



2. APPROACH 



2.1. The Generalized Interpolation Material Point method 

For continuum bodies, the conservation equation for mass is 

^ + pV-v = 0. (1) 

where p is the mass density, v is the velocity, 

In GIMP and MPM, the continuum bodies are discretized with N p material 
particles. Each material particle carries the information of position, velocity, 
temperature, mass, density, Cauchy stress, strain tensor and all other internal state 
variables necessary for the constitutive model. Since the mass of each material particle 
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Figure 1. Contiguous particle GIMP weighting function in one dimension. 



is equal and fixed, Eq.([T]) is automatically satisfied. At each time step, the mass and 
velocities of the material particles are mapped onto the background computational 
mesh(grid). The mapped nodal velocity Vj is obtained through the following equation, 

Yl mi ^i = Yl m p w p N * ( x p) ( 2 ) 
i p 

where m p , v p and x p are the mass, velocity and position of particle p, respectively. N{ 
is the shape function, i and j indexes of node. 

In the early version of MPM, the grid shape function iVj is not smoothed in the 
construction of the weighting function which causes the numerical noise as the material 
points cross computational grid boundaries. Bardenhangen et al [26] presented a family 
of methods named the Generalized Interpolation Material Point (GIMP) methods in 
which the interpolation function are in C 1 (as opposed to MPM, for which are in C ). 
In this paper the following smoothed shape function in C 1 is used, iVj = Q(r x )Q(r y ), 
where r x = \x p — Xi\/L, r y = \y p — yi\/ L, L is the length of cell, $(r) is given in flowing, 
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See also Fig{TJ Eq.Q is a weighting function with support in adjacent cells and the next 
nearest neighbor cells. This specialization has the advantage that it develops weighting 
function in C 1 with a minimal amount of additional complexity. The increased support 
does result in an increase in computational effort that is different from the meshless 
methods which invest a mount of effort to the research of the influence nodes to construct 
the shape function. 

In Eq.(j2]), the consistent mass matrix, rriij, is 

J^m p iV i (xp)JV i (x i ,). (4) 



rriij 
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In practice, we generally replace m^- with a lumped, diagonal mass matrix so that Eq. 
(J2J) becomes 

mjVi = ^JmpVpiVi(xp) (5) 

v 

where lumped mass is 

m i = ^2m p N i (x p ) (6) 

p 

The conversation equation for momentum reads, 
dv „ 

P dt = P ' ( ) 

where cr is the stress tensor and b is the body force. The weak form of Eq.([7]) based on 
the standard procedure used in the FEM [28, 29] can be written 

/ p5v ■ ^dfi + f 5v ■ (a ■ n - t)dr + [ p5v- hdil = (8) 
Jn dt JT t Jn 

where n and t is the outward normal unit and traction vectors on the boundary 

Since the continuum bodies is described with the use of a finite set of material 
particles, the mass density can be written as, 

N p 

p(x) = ^M/(x-x p ) (9) 
p=i 

where 5 is the Dirac delta function with dimension of the inverse of volume. The 
substitution of Eq. (jUJ) into Eq. (jSJ) converts the integral to the sums of quantities 
evaluated at the material particles, namely, 

= (f,) int + (f,) ext (10) 
where mi is the lumped mass and the internal force vector is given by 

Np 

{f i ) iBb = -Y t M p (T'{VN i )/p p (11) 

p 

and the external force vector reads 

N p 

= J^JVibp + tf (12) 
p=i 

where the vector f? is the contact force which is the external nodal force not including 
the body force and is illustrated in the following section. 

A explicit time integrator is used to solve Eq. (flOl for the nodal accelerations, with 
the time step satisfying the stability condition. The critical time step is the ratio of 
the smallest cell size to the wave speed. After the equations of motion are solved on 
the cell nodes, the new nodal values of acceleration are used to update the velocity of 
the material particles. The strain increment for each material particle is determined 
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Figure 2. (See attached Fig2.jpg)The scheme of contact between body p and body q. 

with the use of gradient of the nodal basis function evaluated at the material particle 
position. The corresponding stress increment can be found from the constitutive model. 
The internal state variables can also be completely updated. The computational mesh 
may be discarded, and a new mesh is defined, if desired for the next time step. As a 
result, an effective computational mesh could be chosen for convenience [31]. 

2.2. Contact algorithm 

The GIMP method provides a natural no-slip contact algorithm based on a common 
background mesh. But natural contact algorithm has two drawbacks: firstly, the 
premature contact occurs because the velocities of two bodies are mapped on the same 
nodes though the distance between the bodies may be still two or even more times of the 
length of cell, which may cause the numerical noise of stress; Secondly, it is impossible 
to separate the contacting bodies. Bardenhagen et al have proposed a contact algorithm 
to simulate the interactions of the grains of granular material [40], in which the contact 
between bodies is handled when the velocity field of individual particles in contact differs 
from the single, center-of-mass velocity field in the cell containing contacting particles. 
A multi-mesh mapping scheme is proposed by Hu and Chen [43]. In the multi-mesh 
mapping scheme, each material lies in an individual background mesh rather than in 
the common background one. The meshing process of spur gears is simulated by Hu 
and Chen with their contact algorithm. To avoid interpenetration and allow separation 
in the gear meshing process, the normal velocity of any particle at the contact surface 
is calculated in the common background mesh, while the tangential velocity is found 
based on the corresponding information in respective individual mesh. 

In Bardenhagen's and Chen's contact algorithms the premature contact may still 
occur. In this paper, the criterion of contact consults the distance of the bodies. In the 
improved contact algorithm, if the velocities of body p and body q are mapped on the 
same node i (Seen in FigfS]), the distance between body p and q is calculated. We note 
the distance between body p and q as D^f 3 , which can be calculated as: 

fflM = n pg . r P _ n pg . r ? ( 13) 

where n pq is the normal direction of the contacting interface pointing from body p to 
q. r P and rj are the vectors pointing to node i from body p and body q, respectively. 
n pq ■ r P and n pq ■ r q can be chosen as 

n pg . r P = max j r P . n P? ; m — 1,2, N p and particle p m € body p] (14) 

n pq ■ v\ = min {r^ ■ n pq , m = 1, 2, N q and particle q m G body q] (15) 

where r p and r q are the vectors pointing to node i from particle p m and q m , N p and 
N q are the numbers of particles belonging to body p and q, respectively. 
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The the criterion of contact can be written as 

*r < \, (is) 

where L is the length of cell. If Eq. (fl6l) is satisfied, the velocities of p and q are adjusted 
to new values so that the normal components of them are equal, then the strain and 
stress of all the particles are updated. So, Eq. ffToT) plays the role of a switch function. 
Once bodies p and q contact, they move together along the normal until they separate, 
so the acceleration along the normal of body p is equal to that of q during the course of 
the contact. That is 

a? • nf* = a? • nf (17) 

where af and are the accelerations of bodies p and q at node i, respectively. They 
can be obtained from the Newtonian second law. Also the normal contact force /, 
can be derived. Note that the normal contact force must be nonnegative. So, once /, 
is negative, body p and q is not in contact in the next time step. 

If without friction, the contact algorithm has been finished up to now. In the case 
with friction the Coulomb friction is applied to calculate the tangential contact force. 



nor 

% 

nor 

i 



2. 3. Constitutive model and the Equation of State 

In this paper, the material is modeled using von Mises plasticity with linear harding [3]. 
A plastic model dictates a linear elastic response until a yield criterion is reached. The 
von Mises yield criterion is 3 J 2 — Y 2 = 0, where Y is the plastic yield stress, J 2 is the 
second invariant of s, J 2 = |s : s, s is the deviatoric stress tensor, s = a — Tr[c]/3. 
The linear hardening means that Y increases linearly with the second invariant of the 
plastic strain tensor. If 3J 2 > Y 2 , the increment of equivalent plastic strain de p can 
be calculated as de p = (\f3J 2 — Y)/(3G + E ttm ), where G and E tan are the shear and 
harding modulus, respectively. The increment of the plastic energy can be calculated 
as dW p = de p ■ Y. It is totally translated into the internal energy. 

The pressure P is calculated by using the Mie-Griineissen state of equation which 
can be written as 

P-P H = - E H {V H )\ (18) 

This description consults the Rankine-Hugoniot curve. In Eq. tflBl . Pjj, Vh and Eh are 
pressure, specific volume and energy on the Rankine-Hugoniot curve, respectively. The 
relation between Ph and Vh can be estimated by experiment and can be written as 
f pocgq-^-) 

P H = { (a-d^x^-i)*' VH ~ Vo (19) 
I PlcK^ - I), V H > Vo 

where cq is the sound speed, A is the parameter in the linear ration, U s = cq + XU P , 
U s and Up are shock velocity and velocity of particle, respectively. In this paper, the 
coefficient of Griineissen j(Vh) is taken as a constant and the transformation of specific 
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internal energy E — Eh{Vh) is taken as the plastic energy. Both the shock compression 
and the plastic work cause the increasing of temperature. The increasing of temperature 
from shock compression can be calculated as: 

^Th = cl ■ X(V - V H f _ 7 (V) 

dVff c v [(X- l)Vb- W H f V H H ' 1 ] 

where c v is the specific heat. Eq. (1201) can be resulted with thermal equation and the 
Mie-Griineissen state of equation [45]. The increasing of temperature from plastic work 
can be calculated as: 

dT p = ™t (21) 

c v 

Both the Eq. fl20p and the Eq. (121 j) can be written as the form of increment. In the 
present work the thermal dissipation is not taken into account. Such a treatment is 
reasonable for cases where the propagation speed of shock is much faster than that of 
the thermal dissipation. 

The material constants are chosen to model the energetic crystal high melting 
explosive. The elastic modulo is 11.87GPa, the poisson ratio is 0.25, the density is 
1.9 x 10~ 3 g/mm 3 , the initial yield strength is lOOMPa, the hardening modulus is OMPa, 
the heat capacity is 1250J/(Kg • K). The coefficient of Griineissen is taken as 1.1, A is 
2.6 and the sound speed c is 2740m/s [1]. 



3. Results 



In this paper we focus on the two-dimensional case. The shock wave to material target is 
loaded via the impacting by a second material block with symmetric configuration and 
opposite velocity. The initial shock is along the vertical direction. We apply periodic 
boundary conditions to simulated system in the horizontal direction. As the first step, 
we set a single cavity in the HMX block. Such a simulation model corresponds also to 
a very wide system with a row of cavities being parallel to the impacting plane. We 
study various cases with different cavity radius and with different strengths of shock. 
The influence of existing cavity to the contact of bodies is investigated and a secondary 
impacting is observed. In the cases of strong shock, the influence of cavity size and the 
impact speed to the "hot-spot" is systematically investigated. 



3.1. Validation of the improved contact algorithm 

The newly proposed contact algorithm is validated by simulating the impacting of two 
identical solid bodies having opposite velocities in vertical direction. The initial speed 
of each body is set to be 300m/s. Initially, the distance between the two bodies is set as 
4mm. The width and height of blocks are set to be 40mm and 50mm, respectively. Figj3] 
shows the pressure along the vertical direction calculated by different contact algorithms, 
where the impacting interface is at y — 0. From the numerical results we can see that 
both the GIMP with Bardenhagen's contact algorithm and the MPM without contact 
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Figure 3. (Color online) Pressure distribution along the vertical direction at time 
t = 0.014ms. The impacting interface is at y = 0. Previous contact algorithms present 
unreasonable too high pressure in the vicinity of the impact interface. 



algorithm give an unreasonable and too high pressure in the vicinity of the interface 
corresponding to a severe numerical premature contact (See the dash-dotted and dotted 
lines in the figure). The result by the improved contact algorithm is much better (See 
the red line in the figure). It is well-known that spurious "wall- heating" phenomenon 
is generally difficult to eliminate for computational fluid dynamics and computational 
solid dynamics. Here, our new contact scheme decreases the "wall-heating" phenomenon 
to an extent that is nearly undistinguishable. 

3.2. Dynamical response to weak shock 

Due to the ability of withstanding the deformation of shear and tension, cavities in the 
solid material will not collapse when the shock is very weak. Fig. H] shows a series of 
snapshots for a case with an initial impact speed of 150m/s. The size of blocks is the 
same as section 13.11 and the radius of each cavity is 10 mm. The global procedure can 
be described as follows: 

(a) , The two bodies get contact and a shock wave is loaded to each. The shock 
waves in the two bodies move forwards to the two cavities. The pressure of shock is 
about 900MPa(See Fig^a)). 

(b) , The shocks arrive at the cavities, then rarefactive waves are reflected back to 
material in between the two cavities. The cavities begin to shrink in a nearly isotropic 
fashion(See in Figfjjb)). 

(c) , The rarefactive waves reach the surface of contacting, then the two bodies begin 
to separate(See in Figfjjc)). Note that if there is no cavity inside, the two bodies will 
not separate so early. 

(d) , The two bodies continue to separate until the distance between the two bodies 
reaches the maximal value(See in FigH(d)), then the two bodies begin to move closer 
again because the particles, being above the upper cavity and down the lower cavity, are 
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Figure 4. (See attached Fig4.jpg)Snapshots for the impact of two blocks with a single 
cavity in each. Impact velocity of each block is v — 150m/s. From black to white 
the gray level in the figure shows the increase of local pressure. The unit of pressure 
is MPa. (a) t = 0.004ms, (b) t = 0.008ms, (c) t = 0.014ms, (d) t = 0.025ms, (e) 
t = 0.06ms, (f) t = 0.18ms. 
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Figure 5. The separation distance versus time t for different impacting speeds shown 
in the legend. D represents half of the distance between the two impacting surfaces. 
Before the occurrence of the secondary impact, D has a maximum value D max which 
increases with the impacting speed. 

still moving towards the contacting surface, which draws the separating bodies back. 

(e) , The approaching two bodies collide again. They stay together until the 
rarefactive waves reflected from the lowermost and the uppermost free surfaces arrive 
at the contacting surface. Then, they will separate again (See in FigJUe)). 

(f) , After e, the two bodies separate. This is the final separation. The terminal 
configuration of cavities are still approximately circular. From the whole procedure, we 
can see that, due to the existence of cavities, the two bodies get a secondary impacting. 

In order to comprehend more clearly the phenomena of the secondary impacting, we 
alter the impacting speed(in the limit of no collapse) and the size of cavity. Figj5] shows 
the effects of the initial impacting speed to the maximum separation distance between 
the two bodies after the first impacting. In Figj5l D represents half of the distance 
between the two impacting surfaces. It varies with time t. The initial impacting speed 
is set as lOOm/s, 150m/s and 200m/s, respectively. In all cases the phenomena of the 
secondary impacting are observed. The maximum value of D is denoted by -D max . The 
larger the initial impacting speed, the larger the value of D max . 

Figj6] shows the effects of cavity size on the separation distance between the two 
impacting bodies. Here D and D max have the same meanings as in Figj5l The initial 
radius of each cavity is set as 6mm, 8mm and 10mm, respectively. It can be seen that 
the larger the radius of cavity, the larger the value of -D max . That is to say, the larger the 
cavities, the easier to observe the phenomenon of the secondary impacting. Naturally, 
if the radius of cavity diminishes to zero, the secondary impact disappears. 
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Figure 6. The distance of separation versus time t for different cavity sizes shown in 
the legend. Here D and -D max have the same meanings as in FigEl 

Figure 7. (See attached Fig7.jpg) Snapshots for the impacting process of two blocks 
with a single cavity. The initial radius of cavity 10mm, impact speed v = 300m/s. From 
blue to red, the local pressure increases. The unit of pressure is Mpa. t = 0.008ms in 
(a), t = 0.012ms in (b), t = 0.02ms in (c), t = 0.024ms in (d). 



3.3. Dynamical response to strong shock 

If the shock is strong, asymmetric collapse of the cavity will occur. When the radius 
of cavity is set to be 6mm, 8mm and 10mm,respectively, the corresponding minimum 
impact speed v for the cavity to collapse asymmetrically is 120m/s, 150m/s and 200m/s, 
respectively That is to say, the smaller the radius of cavity, the easier for the cavity 
to collapse asymmetrically. This is because larger curvature is easier to accumulate 
particles. The course of the collapse of the upper cavity can be described as follows: 
When the shock reaches the cavity, the upstream side of cavity will reflect a rarefactive 
wave which propagates downwards. See FigJT^a). This causes the pressure in the region 
close to the upstream side of the cavity is lower than around. Then, the neighboring 
particles accumulate towards to the lower side of cavity, which will accelerate the upward 
speed of particles there. There are two factors that influence the deformation of the 
cavity. One is the the accumulation of the particles to the lower side of the cavity, 
the other is the resistance of the solid material to shear. The competition of the two 
factors determines how the cavity deforms. Obviously, as the impact speed increases, 
the resulted shock becomes stronger, then the influence of the accumulation of the 
particles becomes stronger. If the strength of the material close to the lower side of the 
cavity cannot withstand the accumulation of particles, the upward speed of particles 
will become larger and larger, and the cavity will collapse asymmetrically. The upward 
speed of these particles in FigJT^a) is about 550m/s, increases to 844m/s in FigfT^b), 
reaches lOOOm/s and a jet can be observed in FigJT^c), after that, these particles rush 
to the upper side of cavity at a speed about lOOOm/s. At the time of FigJT^d) the jet 
just impacts the upper side of cavity, which causes of a distinct increase of pressure. 
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Figure 8. The maximum relative velocity versus the initial impacting speed. The 
radius of cavity is set as 6mm, 8mm and 10mm, respectively. 




0.35 



Figure 9. The maximum relative velocity versus the radius of cavity. The impact 
speed is set as 500m/s. 



The course of collapsing of the lower cavity is similar. 

To quantitatively describe the collapsing procedure of the shocked cavity, we 
measure the maximum speed of particles being close to the upper-steam wall relative 
to that being close to the down-stream wall. FigJH] shows the maximum relative 
velocity versus to the impacting speed. For the investigated cases, the former increases 
nonlinearly with the impacting speed and approaches to a constant; the smaller the 
initial radius, the larger the maximum relative velocity. A specific case with an 
impacting speed v = 500 m/s is shown in FiglHl The approached constant is about 
3.58Km/s. 

Also, the phenomena of the secondary impacting are observed in cases the with 
strong shock. Fig{10] shows the effects of initial impacting speed v on the separation 
distance D between the two impacting bodies. The initial impacting speed v is 
set as 230m/s, 300m/s, 500m/s and 700m/s, respectively. Contrast to the cases 
with weak shock, -D max doesn't increase monotonously with the initial impact speed. 
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Figure 10. (Color online) The separation distance versus time t for different impacting 
speeds with strong shock. Hear D and has the same meaning as FigE] 



Figure 11. (See attached Figff.jpg) Snapshots for the impact of two blocks with a 
single cavity in each. Impact velocity v — 300m/s. From blue to red, the color in 
the figure shows the increase of local temperature. The unit of temperature is K. (a) 
t = 0.004ms, (b) t = 0.006ms, (c) t = 0.014ms, (d) t = 0.020ms, (e) t = 0.025ms, (f) 
t = 0.030ms. 



When the initial impact speed is not very large (lower than 300m/s), -D max increases 
correspondingly; but when the impact speed is large (higher than 300m/s), -D max 
decreases. The critical case with v = 300 m/s is referred to FigJTl A physical explanation 
based on the case with R = 10 mm is as below. When the impacting speed v is lower 
than 300 m/s (see Fig. @J, the cavities shrink but without jetting. The strength of the 
raref active wave increases with the initial impacting speed v. When v is larger than 
300 m/s, the cavities collapse and cave-in which dissipates a considerable amount of the 
kinetic energy. The amount of the dissipated energy exceeds the increasing of the total 
kinetic energy via increasing v. Therefore, the strength of the rarefactive waves do not 
increase with v but decrease. 

34. On "hot-spots" 

The reactivity of porous energetic material depends greatly on the nature of the "hot- 
spots" formed by shocks as they move through the material. As observed in the gas-gun 
experiments [44], coarse HMX produce "hot-spots" that are large enough to persist for 
a long time whereas fine material seems to produce "hot-spots" that are smaller in size 
and dissipate quickly. The threshold-to-initiation is the limit where exothermic chemical 
energy release is balanced by energy dissipated away from the "hot-spot" reaction. In 
the work of Baer, et al [1], the mesoscale processes of consolidation, deformation and 
reaction of shocked porous energetic materials are studied using shock physics analysis 
of impact on a collection of discrete HMX "crystals". In our work, we focus on the 
"hot-spots" produced in the shocked HMX material with cavities. 
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Figure 12. The impact speed dependence of the temperature of the "hot spot". The 
radius of cavities is set as 6mm, 8mm and 10mm, respectively. The unit of temperature 
is K. 

FigJTTI shows the global procedure of the production of the "hot-spots". The impact 
speed is 300m/s, the initial sizes of blocks are the same as those in section I3TT] and the 
initial radius of every cavity is 10mm. The initial temperature of the simulated system 
is 300K. In order to exhibit clearly, only the upper block is shown. The procedure can 
be described as follows, 

(a) , The two blocks collide and the shocks begin to propagate towards the cavities. 
The temperature in the shocked pure HMX is about 310K(See in FigJTlTa)). 

(b) , The shocks reach the cavities, then rarefactive waves are reflected back from 
the upstream boundaries of the cavities, the temperatures in the regions which are first 
shocked and then rarefacted by waves decrease a little (See in Fig JTiT b)). 

(c) , The shocks continue to propagate. The sizes of the cavities reduce as time goes 
on, the vertical deformation is more serious than the horizontal one. Opposite to the 
last item, the temperatures in the regions neighboring to the upstream sides of cavities 
increase. The reason is that the plastic work increases continuously (See in Fig flTT c)). 

(d) , A jet appears in each cavity, and the temperature of jet material increases as 
time goes on(See in Fig JTlT d)). 

(e) , The jet material impacts the downstream side of the cavity, which cause a 
distinct increase of temperature to about 520K (See in Fig JTlT e)). 

(f) , The cavities continues to collapse until they are crammed, a "hot-spot" with 
the temperature about 520K is produced in each colliding region (See in FigfTlTf)). 

Fig{12] shows the temperatures of "hot-spot" produced by different impact speeds 
for different sizes of cavities. It can be seen that the temperature increases nearly 
parabolically with the impacting speed v for a fixed size of cavity, which is consistent 
with the Hugoniot relation, P = po(co + Ait), where u = v/2 is the particle speed 
after shock. At the same time, plastic work increases also nearly parabolically with the 
impacting speed v. Also, it can be found that the temperature of "hot-spot" produced 
in smaller cavity is higher than that in larger ones with the same impact speed. As 
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analyzed in section 13.31 the maximum relative speed produced by smaller cavities is 
larger than that by larger ones, which cause higher temperature, for the investigated 
cases. 

4. Conclusion 

The dynamical response of high melting explosive with cavities under shock is 
investigated by the generalized interpolation material point method where the criterion 
for contacting is improved. An elastic-to-plastic and thermal-dynamical model with 
the Mie-Griineissen equation of state consulting the Rankine-Hugoniot curve is used to 
simulate the mechanical and thermal behaviors. A phenomena of secondary impacting 
is observed for the case of colliding of two bodies with a single cavity in each. For 
weak shocks, the cavities shrink in a nearly symmetric way, the phenomena of secondary 
impacting becomes more distinctly as the impact speed increases. But for strong shocks, 
the situation becomes complex because the asymmetric collapse of cavity influences the 
reflection of shock. The course for the collapse of cavity is studied for various cavity 
sizes. For the checked cases, smaller cavities collapse more easily. Our numerical 
results show that the existence of cavities greatly help the creation of "hot-spots" 
with a much higher temperature, which is very important for igniting of the explosive 
materials. The influence of various factors, including impact speed, size of cavity, to the 
temperature of "hot-spots" is investigated. The temperature of "hot-spots" increases 
nearly parabolically with the impact speed for a fixed cavity size. The smaller the 
cavity, the higher the temperature of "hot-spots", for cases investigated. If the cavity 
is too small, the effects of cavity collapse becomes very weak. Future work includes the 
incorporation of thermal diffusion, strain rate effects, etc into the simulations, and cases 
with cavities filled by gases. 
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